XBhat <- as.matrix(newX)%*%Bhat
prob <- 1/(1+exp(-XBhat))
decay.plot <- qplot(x=decayseq, y=prob, geom="line", xlab="t", ylab="Pr(Alliance Gain | t, median(X)")
decay.plot <- decay.plot + ggtitle("Probability of Gaining a New Alliance\nAfter an Alliance Violation")+theme_bw()+theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank())
ggsave(decay.plot, file="Figure_A2.pdf", width= 6, height=6)
# clear memory in R
rm(list=ls())
# load libraries
library(ggplot2)
library(dplyr)
library(arm)
library(foreign)
library(games)
library(sandwich)
library(lmtest)
library(MASS)
install.packages("ggplot2")
# clear memory in R
rm(list=ls())
# load libraries
library(ggplot2)
library(dplyr)
library(arm)
library(foreign)
library(games)
library(sandwich)
library(lmtest)
library(MASS)
install.packages("dplyr")
# clear memory in R
rm(list=ls())
# load libraries
library(ggplot2)
library(dplyr)
library(arm)
library(foreign)
library(games)
library(sandwich)
library(lmtest)
library(MASS)
#read in data
alliance <- read.dta("alliance_replication.dta")
############################################
# Create Figure 1 from the main manuscript #
############################################
fig1 <- ggplot(data= subset(alliance, regyrvio_avg > 0), aes(x=year, y=regyrvio_avg))
fig1 <- fig1 + geom_point() + geom_smooth()
fig1 <- fig1 + scale_x_continuous(breaks=seq(1920, 2000, by=20))
fig1 <- fig1 + facet_grid( . ~ region1)
fig1 <- fig1 + theme_bw() + xlab("Year") + ylab("Percent of Countries Ending \n an Alliance Through Violation")
ggsave(fig1, file="Figure_1.pdf", width=10, height = 3)
############################################
# Create Figure 2 from the main manuscript #
############################################
# function to calculate clustered standard errors
# same as the "cluster"" command in stata
clx <- function(fm, dfcw, cluster){
M <- length(unique(cluster))
N <- length(cluster)
dfc <- (M/(M-1))*((N-1)/(N-fm$rank))
u <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum))
vcovCL <- dfc*sandwich(fm, meat=crossprod(u)/N)*dfcw
coeftest(fm, vcovCL)
}
# clx2 returns a variance covariance matrix for cluster robust se's
clx2 <- function(fm, dfcw, cluster){
M <- length(unique(cluster))
N <- length(cluster)
dfc <- (M/(M-1))*((N-1)/(N-fm$rank))
u <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum))
vcovCL <- dfc*sandwich(fm, meat=crossprod(u)/N)*dfcw
vcovCL
}
# Fit the 6 model's from Figure 2 in the main manuscript (from left to right)
# The coeffiecients and standard errors from these models will be used to create Figure 2.
# The standardize command rescales the regression inputs so that they
# are centered at 0 and normalized by 2 sd's
# The coeffients for these normalized inputs form the basis for our
# graph "Figure 2" in the main manuscript.
# Create Matrix to store coefficients
coef.matrix <- matrix(nrow = 6, ncol = 6)
# Create Matrix to store clustered standar errors
se.matrix <- matrix(nrow = 6, ncol = 6)
# Model 1 and 2 comparing our measures with state clustered SE
# but no time dynamics
#Model 1
syrfe.cumvio.notime <- standardize(glm(alliance_gain ~  cum_vio_dum + majpow + democ_dum + lag_num + I(lag_num^2) + as.factor(region1), family= binomial(link=logit),  data=alliance, x=TRUE))
#estimated coefficients and se's using cluster robust se's assuming correlation
# within state observations
clustse.syrfe.cumvio.notime <- clx(syrfe.cumvio.notime, 1, alliance$state[!is.na(alliance$lag_num)])
# store model coeffients and clustered se's
coef.matrix[2:6,1] <- clustse.syrfe.cumvio.notime[2:6,1]
se.matrix[2:6,1] <- clustse.syrfe.cumvio.notime[2:6,2]
#variance covariance matrix
vcov.clustse.syrfe.cumvio.notime <- clx2(syrfe.cumvio.notime, 1, alliance$state[!is.na(alliance$lag_num)])
#Model 2
syrfe.rep.notime <- standardize(glm(alliance_gain ~  rep_dum_mem_indef + majpow + democ_dum + lag_num + I(lag_num^2) + as.factor(region1), family= binomial(link=logit),  data=alliance, x=TRUE))
clustse.syrfe.rep.notime <- clx(syrfe.rep.notime, 1, alliance$state[!is.na(alliance$lag_num)])
# store model coeffients and clustered se's
coef.matrix[c(1,3:6),4] <- clustse.syrfe.rep.notime[2:6,1]
se.matrix[c(1,3:6),4] <- clustse.syrfe.rep.notime[2:6,2]
### Perform Clarke test for
### non-nested models on Model's 1 and 2 as described on
### page 29 of the manuscript
clarke(syrfe.cumvio.notime, syrfe.rep.notime)
# Model 3 and 4 comparing our measures with state clustered SE
# plus time dynamics
#Model 3
syrfe.cumvio.time <- standardize(glm(alliance_gain ~  cum_vio_dum+ majpow + democ_dum + lag_num  + I(lag_num^2)  + t + I(t^2) + I(t^3)  + as.factor(region1), family= binomial(link=logit),  data=alliance))
clustse.syrfe.cumvio.time <- clx(syrfe.cumvio.time, 1, alliance$state[!is.na(alliance$lag_num)])
# store model coeffients and clustered se's
coef.matrix[2:6,2] <- clustse.syrfe.cumvio.time[2:6,1]
se.matrix[2:6,2] <- clustse.syrfe.cumvio.time[2:6,2]
#Model 4
syrfe.rep.time <- standardize(glm(alliance_gain ~  rep_dum_mem_indef + majpow + democ_dum + lag_num  + I(lag_num^2)  + t + I(t^2) + I(t^3)  + as.factor(region1), family= binomial(link=logit),  data=alliance))
clustse.syrfe.rep.time <- clx(syrfe.rep.time, 1, alliance$state[!is.na(alliance$lag_num)])
# store model coeffients and clustered se's
coef.matrix[c(1,3:6),5] <- clustse.syrfe.rep.time[2:6,1]
se.matrix[c(1,3:6),5] <- clustse.syrfe.rep.time[2:6,2]
#Model 5 and 6  add a decay term
# Model 5
syrfe.cumvio.decay <- glm(alliance_gain ~  cum_vio_dum + majpow + democ_dum + lag_num  + I(lag_num^2)  + t + I(t^2) + I(t^3)  + decay_gen + I(decay_gen^2) + I(decay_gen^3)+ as.factor(region1), family= binomial(link=logit),  data=alliance)
clustse.syrfe.cumvio.decay <- clx(syrfe.cumvio.decay, 1, alliance$state[!is.na(alliance$lag_num)])
# store model coeffients and clustered se's
coef.matrix[2:6,3] <- clustse.syrfe.cumvio.decay[2:6,1]
se.matrix[2:6,3] <- clustse.syrfe.cumvio.decay[2:6,2]
# Model 6
# unnormalized model for later
syrfe.rep.decay <- glm(alliance_gain ~  rep_dum_mem_indef + majpow + democ_dum + lag_num  + I(lag_num^2)  + t + I(t^2) + I(t^3)  + decay_gen + I(decay_gen^2) + I(decay_gen^3)+ as.factor(region1), family= binomial(link=logit),  data=alliance)
clustse.syrfe.rep.decay <- clx(syrfe.rep.decay, 1, alliance$state[!is.na(alliance$lag_num)])
# normalized for Figure2
syrfe.rep.decay.norm <- standardize(glm(alliance_gain ~  rep_dum_mem_indef + majpow + democ_dum + lag_num  + I(lag_num^2)  + t + I(t^2) + I(t^3)  + decay_gen + I(decay_gen^2) + I(decay_gen^3)+ as.factor(region1), family= binomial(link=logit),  data=alliance))
clustse.syrfe.rep.decay.norm <- clx(syrfe.rep.decay.norm, 1, alliance$state[!is.na(alliance$lag_num)])
# store model coeffients and clustered se's
coef.matrix[c(1,3:6),6] <- clustse.syrfe.rep.decay.norm[2:6,1]
se.matrix[c(1,3:6),6] <- clustse.syrfe.rep.decay.norm[2:6,2]
# Run additional code to make  Figure 2
source("table_graph_code.R")
######################################
# Create Figure 3 from the main text #
######################################
# This graphs the estimated effect of reputation
attach(alliance)
sims <- 5000
Bhat <- as.vector(clustse.syrfe.rep.decay[,1])
par.draws <- mvrnorm(sims, mu=Bhat, Sigma=clx2(syrfe.rep.decay, 1, alliance$state[!is.na(alliance$lag_num)]))
# set all variables to their median
X1 <- cbind(1, 0, median(majpow, na.rm=T), median(democ_dum, na.rm=T), median(lag_num, na.rm=T), median(lag_num, na.rm=T)^2, median(t, na.rm=T), median(t, na.rm=T)^2, median(t, na.rm=T)^3, median(decay_gen, na.rm=T), median(decay_gen, na.rm=T)^2, median(decay_gen, na.rm=T)^3, 0,0,0,0)
est.fx <- c()
for(i in 6:81){
X1[2] <- i
mu <- 1/(1+exp(-(par.draws%*%t(X1))))
mean.fx <- mean(mu)
sd.fx <- sd(mu)
ci.fx <- quantile(mu, c(.025, 0.975))
est.fx <- rbind(est.fx, c(mean.fx,sd.fx,ci.fx))
}
est.fx <- as.data.frame(est.fx)
names(est.fx) <- c("mean","sd","ci.lo", "ci.hi")
fx.plot <- ggplot(data=est.fx, aes(x=seq(6,81,1), y=mean))+
geom_ribbon(aes(x=seq(6,81,1), ymin=ci.lo, ymax = ci.hi, alpha=.3))+
geom_line(size = 2)+
xlab("Unreliability Score")+ylab("Pr(Alliance Gain | 1 Violation)")
fx.plot <- fx.plot + ggtitle("The Probability of Gaining a new Alliance\nas a Function of Unreliability") +theme_bw() +  theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank(), legend.position="none")
ggsave(fx.plot, file="Figure_3.pdf", width= 6, height =6)
####################################
# Create figure A1 in the appendix #
####################################
timeseq <- 0:max(alliance$t, na.rm=T)
# Set all variables to their median
newX2 <- cbind(1, median(rep_dum_mem_indef, na.rm=T), median(majpow, na.rm=T), median(democ_dum, na.rm=T), median(lag_num, na.rm=T), median(lag_num, na.rm=T)^2, timeseq, timeseq^2,timeseq^3,0,0,0, 0,0,0,0)
Bhat2 <- as.vector(clustse.syrfe.rep.decay[,1])
XBhat2 <- as.matrix(newX2)%*%Bhat2
prob2 <- 1/(1+exp(-XBhat2))
time.plot <- qplot(x=timeseq, y=prob2, geom="line", xlab="t", ylab="Pr(Alliance Gain | t, median(X)")
time.plot <- time.plot + ggtitle("Probability of Gaining a New Alliance\nAfter an Alliance Gain") +theme_bw() +theme( panel.grid.minor = element_blank(), panel.grid.major = element_blank())
ggsave(time.plot, file="Figure_A1.pdf", height=6, width=6)
####################################
# Create Figure A2 in the appendix #
####################################
#graph decay
decayseq <- 0:max(alliance$decay_gen, na.rm=T)
#Set all variables to their median
newX <- cbind(1, median(rep_dum_mem_indef, na.rm=T), median(majpow, na.rm=T), median(democ_dum, na.rm=T), median(lag_num, na.rm=T), median(lag_num, na.rm=T)^2,median(t, na.rm=T), median(t, na.rm=T)^2, median(t, na.rm=T)^3, decayseq, decayseq^2,decayseq^3, 0,0,0,0)
Bhat <- as.vector(clustse.syrfe.rep.decay[,1])
XBhat <- as.matrix(newX)%*%Bhat
prob <- 1/(1+exp(-XBhat))
decay.plot <- qplot(x=decayseq, y=prob, geom="line", xlab="t", ylab="Pr(Alliance Gain | t, median(X)")
decay.plot <- decay.plot + ggtitle("Probability of Gaining a New Alliance\nAfter an Alliance Violation")+theme_bw()+theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank())
ggsave(decay.plot, file="Figure_A2.pdf", width= 6, height=6)
# clear memory in R
rm(list=ls())
# load libraries
library(ggplot2)
library(dplyr)
library(arm)
library(foreign)
library(games)
library(sandwich)
library(lmtest)
library(MASS)
#read in data
alliance <- read.dta("alliance_replication.dta")
############################################
# Create Figure 1 from the main manuscript #
############################################
fig1 <- ggplot(data= subset(alliance, regyrvio_avg > 0), aes(x=year, y=regyrvio_avg))
fig1 <- fig1 + geom_point() + geom_smooth()
fig1 <- fig1 + scale_x_continuous(breaks=seq(1920, 2000, by=20))
fig1 <- fig1 + facet_grid( . ~ region1)
fig1 <- fig1 + theme_bw() + xlab("Year") + ylab("Percent of Countries Ending \n an Alliance Through Violation")
ggsave(fig1, file="Figure_1.pdf", width=10, height = 3)
############################################
# Create Figure 2 from the main manuscript #
############################################
# function to calculate clustered standard errors
# same as the "cluster"" command in stata
clx <- function(fm, dfcw, cluster){
M <- length(unique(cluster))
N <- length(cluster)
dfc <- (M/(M-1))*((N-1)/(N-fm$rank))
u <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum))
vcovCL <- dfc*sandwich(fm, meat=crossprod(u)/N)*dfcw
coeftest(fm, vcovCL)
}
# clx2 returns a variance covariance matrix for cluster robust se's
clx2 <- function(fm, dfcw, cluster){
M <- length(unique(cluster))
N <- length(cluster)
dfc <- (M/(M-1))*((N-1)/(N-fm$rank))
u <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum))
vcovCL <- dfc*sandwich(fm, meat=crossprod(u)/N)*dfcw
vcovCL
}
# Fit the 6 model's from Figure 2 in the main manuscript (from left to right)
# The coeffiecients and standard errors from these models will be used to create Figure 2.
# The standardize command rescales the regression inputs so that they
# are centered at 0 and normalized by 2 sd's
# The coeffients for these normalized inputs form the basis for our
# graph "Figure 2" in the main manuscript.
# Create Matrix to store coefficients
coef.matrix <- matrix(nrow = 6, ncol = 6)
# Create Matrix to store clustered standar errors
se.matrix <- matrix(nrow = 6, ncol = 6)
# Model 1 and 2 comparing our measures with state clustered SE
# but no time dynamics
#Model 1
syrfe.cumvio.notime <- standardize(glm(alliance_gain ~  cum_vio_dum + majpow + democ_dum + lag_num + I(lag_num^2) + as.factor(region1), family= binomial(link=logit),  data=alliance, x=TRUE))
#estimated coefficients and se's using cluster robust se's assuming correlation
# within state observations
clustse.syrfe.cumvio.notime <- clx(syrfe.cumvio.notime, 1, alliance$state[!is.na(alliance$lag_num)])
# store model coeffients and clustered se's
coef.matrix[2:6,1] <- clustse.syrfe.cumvio.notime[2:6,1]
se.matrix[2:6,1] <- clustse.syrfe.cumvio.notime[2:6,2]
#variance covariance matrix
vcov.clustse.syrfe.cumvio.notime <- clx2(syrfe.cumvio.notime, 1, alliance$state[!is.na(alliance$lag_num)])
#Model 2
syrfe.rep.notime <- standardize(glm(alliance_gain ~  rep_dum_mem_indef + majpow + democ_dum + lag_num + I(lag_num^2) + as.factor(region1), family= binomial(link=logit),  data=alliance, x=TRUE))
clustse.syrfe.rep.notime <- clx(syrfe.rep.notime, 1, alliance$state[!is.na(alliance$lag_num)])
# store model coeffients and clustered se's
coef.matrix[c(1,3:6),4] <- clustse.syrfe.rep.notime[2:6,1]
se.matrix[c(1,3:6),4] <- clustse.syrfe.rep.notime[2:6,2]
### Perform Clarke test for
### non-nested models on Model's 1 and 2 as described on
### page 29 of the manuscript
clarke(syrfe.cumvio.notime, syrfe.rep.notime)
# Model 3 and 4 comparing our measures with state clustered SE
# plus time dynamics
#Model 3
syrfe.cumvio.time <- standardize(glm(alliance_gain ~  cum_vio_dum+ majpow + democ_dum + lag_num  + I(lag_num^2)  + t + I(t^2) + I(t^3)  + as.factor(region1), family= binomial(link=logit),  data=alliance))
clustse.syrfe.cumvio.time <- clx(syrfe.cumvio.time, 1, alliance$state[!is.na(alliance$lag_num)])
# store model coeffients and clustered se's
coef.matrix[2:6,2] <- clustse.syrfe.cumvio.time[2:6,1]
se.matrix[2:6,2] <- clustse.syrfe.cumvio.time[2:6,2]
#Model 4
syrfe.rep.time <- standardize(glm(alliance_gain ~  rep_dum_mem_indef + majpow + democ_dum + lag_num  + I(lag_num^2)  + t + I(t^2) + I(t^3)  + as.factor(region1), family= binomial(link=logit),  data=alliance))
clustse.syrfe.rep.time <- clx(syrfe.rep.time, 1, alliance$state[!is.na(alliance$lag_num)])
# store model coeffients and clustered se's
coef.matrix[c(1,3:6),5] <- clustse.syrfe.rep.time[2:6,1]
se.matrix[c(1,3:6),5] <- clustse.syrfe.rep.time[2:6,2]
#Model 5 and 6  add a decay term
# Model 5
syrfe.cumvio.decay <- standardize(alliance_gain ~  cum_vio_dum + majpow + democ_dum + lag_num  + I(lag_num^2)  + t + I(t^2) + I(t^3)  + decay_gen + I(decay_gen^2) + I(decay_gen^3)+ as.factor(region1), family= binomial(link=logit),  data=alliance)
clustse.syrfe.cumvio.decay <- clx(syrfe.cumvio.decay, 1, alliance$state[!is.na(alliance$lag_num)])
# store model coeffients and clustered se's
coef.matrix[2:6,3] <- clustse.syrfe.cumvio.decay[2:6,1]
se.matrix[2:6,3] <- clustse.syrfe.cumvio.decay[2:6,2]
# Model 6
# unnormalized model for later
syrfe.rep.decay <- glm(alliance_gain ~  rep_dum_mem_indef + majpow + democ_dum + lag_num  + I(lag_num^2)  + t + I(t^2) + I(t^3)  + decay_gen + I(decay_gen^2) + I(decay_gen^3)+ as.factor(region1), family= binomial(link=logit),  data=alliance)
clustse.syrfe.rep.decay <- clx(syrfe.rep.decay, 1, alliance$state[!is.na(alliance$lag_num)])
# normalized for Figure2
syrfe.rep.decay.norm <- standardize(glm(alliance_gain ~  rep_dum_mem_indef + majpow + democ_dum + lag_num  + I(lag_num^2)  + t + I(t^2) + I(t^3)  + decay_gen + I(decay_gen^2) + I(decay_gen^3)+ as.factor(region1), family= binomial(link=logit),  data=alliance))
clustse.syrfe.rep.decay.norm <- clx(syrfe.rep.decay.norm, 1, alliance$state[!is.na(alliance$lag_num)])
# store model coeffients and clustered se's
coef.matrix[c(1,3:6),6] <- clustse.syrfe.rep.decay.norm[2:6,1]
se.matrix[c(1,3:6),6] <- clustse.syrfe.rep.decay.norm[2:6,2]
# Run additional code to make  Figure 2
source("table_graph_code.R")
######################################
# Create Figure 3 from the main text #
######################################
# This graphs the estimated effect of reputation
attach(alliance)
sims <- 5000
Bhat <- as.vector(clustse.syrfe.rep.decay[,1])
par.draws <- mvrnorm(sims, mu=Bhat, Sigma=clx2(syrfe.rep.decay, 1, alliance$state[!is.na(alliance$lag_num)]))
# set all variables to their median
X1 <- cbind(1, 0, median(majpow, na.rm=T), median(democ_dum, na.rm=T), median(lag_num, na.rm=T), median(lag_num, na.rm=T)^2, median(t, na.rm=T), median(t, na.rm=T)^2, median(t, na.rm=T)^3, median(decay_gen, na.rm=T), median(decay_gen, na.rm=T)^2, median(decay_gen, na.rm=T)^3, 0,0,0,0)
est.fx <- c()
for(i in 6:81){
X1[2] <- i
mu <- 1/(1+exp(-(par.draws%*%t(X1))))
mean.fx <- mean(mu)
sd.fx <- sd(mu)
ci.fx <- quantile(mu, c(.025, 0.975))
est.fx <- rbind(est.fx, c(mean.fx,sd.fx,ci.fx))
}
est.fx <- as.data.frame(est.fx)
names(est.fx) <- c("mean","sd","ci.lo", "ci.hi")
fx.plot <- ggplot(data=est.fx, aes(x=seq(6,81,1), y=mean))+
geom_ribbon(aes(x=seq(6,81,1), ymin=ci.lo, ymax = ci.hi, alpha=.3))+
geom_line(size = 2)+
xlab("Unreliability Score")+ylab("Pr(Alliance Gain | 1 Violation)")
fx.plot <- fx.plot + ggtitle("The Probability of Gaining a new Alliance\nas a Function of Unreliability") +theme_bw() +  theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank(), legend.position="none")
ggsave(fx.plot, file="Figure_3.pdf", width= 6, height =6)
####################################
# Create figure A1 in the appendix #
####################################
timeseq <- 0:max(alliance$t, na.rm=T)
# Set all variables to their median
newX2 <- cbind(1, median(rep_dum_mem_indef, na.rm=T), median(majpow, na.rm=T), median(democ_dum, na.rm=T), median(lag_num, na.rm=T), median(lag_num, na.rm=T)^2, timeseq, timeseq^2,timeseq^3,0,0,0, 0,0,0,0)
Bhat2 <- as.vector(clustse.syrfe.rep.decay[,1])
XBhat2 <- as.matrix(newX2)%*%Bhat2
prob2 <- 1/(1+exp(-XBhat2))
time.plot <- qplot(x=timeseq, y=prob2, geom="line", xlab="t", ylab="Pr(Alliance Gain | t, median(X)")
time.plot <- time.plot + ggtitle("Probability of Gaining a New Alliance\nAfter an Alliance Gain") +theme_bw() +theme( panel.grid.minor = element_blank(), panel.grid.major = element_blank())
ggsave(time.plot, file="Figure_A1.pdf", height=6, width=6)
####################################
# Create Figure A2 in the appendix #
####################################
#graph decay
decayseq <- 0:max(alliance$decay_gen, na.rm=T)
#Set all variables to their median
newX <- cbind(1, median(rep_dum_mem_indef, na.rm=T), median(majpow, na.rm=T), median(democ_dum, na.rm=T), median(lag_num, na.rm=T), median(lag_num, na.rm=T)^2,median(t, na.rm=T), median(t, na.rm=T)^2, median(t, na.rm=T)^3, decayseq, decayseq^2,decayseq^3, 0,0,0,0)
Bhat <- as.vector(clustse.syrfe.rep.decay[,1])
XBhat <- as.matrix(newX)%*%Bhat
prob <- 1/(1+exp(-XBhat))
decay.plot <- qplot(x=decayseq, y=prob, geom="line", xlab="t", ylab="Pr(Alliance Gain | t, median(X)")
decay.plot <- decay.plot + ggtitle("Probability of Gaining a New Alliance\nAfter an Alliance Violation")+theme_bw()+theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank())
ggsave(decay.plot, file="Figure_A2.pdf", width= 6, height=6)
rm(list=ls())
setwd("C:\\Users\\hanse\\Dropbox\\Journal (II)\\Replication Files\\Narang\\Final_Replication_Files\\Final_Replication_Files\\")
# load libraries
library(ggplot2)
library(dplyr)
library(arm)
library(foreign)
library(games)
library(sandwich)
library(lmtest)
library(MASS)
#read in data
alliance <- read.dta("alliance_replication.dta")
############################################
# Create Figure 1 from the main manuscript #
############################################
fig1 <- ggplot(data= subset(alliance, regyrvio_avg > 0), aes(x=year, y=regyrvio_avg))
fig1 <- fig1 + geom_point() + geom_smooth()
fig1 <- fig1 + scale_x_continuous(breaks=seq(1920, 2000, by=20))
fig1 <- fig1 + facet_grid( . ~ region1)
fig1 <- fig1 + theme_bw() + xlab("Year") + ylab("Percent of Countries Ending \n an Alliance Through Violation")
ggsave(fig1, file="Figure_1.pdf", width=10, height = 3)
############################################
# Create Figure 2 from the main manuscript #
############################################
# function to calculate clustered standard errors
# same as the "cluster"" command in stata
clx <- function(fm, dfcw, cluster){
M <- length(unique(cluster))
N <- length(cluster)
dfc <- (M/(M-1))*((N-1)/(N-fm$rank))
u <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum))
vcovCL <- dfc*sandwich(fm, meat=crossprod(u)/N)*dfcw
coeftest(fm, vcovCL)
}
# clx2 returns a variance covariance matrix for cluster robust se's
clx2 <- function(fm, dfcw, cluster){
M <- length(unique(cluster))
N <- length(cluster)
dfc <- (M/(M-1))*((N-1)/(N-fm$rank))
u <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum))
vcovCL <- dfc*sandwich(fm, meat=crossprod(u)/N)*dfcw
vcovCL
}
# Fit the 6 model's from Figure 2 in the main manuscript (from left to right)
# The coeffiecients and standard errors from these models will be used to create Figure 2.
# The standardize command rescales the regression inputs so that they
# are centered at 0 and normalized by 2 sd's
# The coeffients for these normalized inputs form the basis for our
# graph "Figure 2" in the main manuscript.
# Create Matrix to store coefficients
coef.matrix <- matrix(nrow = 6, ncol = 6)
# Create Matrix to store clustered standar errors
se.matrix <- matrix(nrow = 6, ncol = 6)
# Model 1 and 2 comparing our measures with state clustered SE
# but no time dynamics
#Model 1
syrfe.cumvio.notime <- standardize(glm(alliance_gain ~  cum_vio_dum + majpow + democ_dum + lag_num + I(lag_num^2) + as.factor(region1), family= binomial(link=logit),  data=alliance, x=TRUE))
#estimated coefficients and se's using cluster robust se's assuming correlation
# within state observations
clustse.syrfe.cumvio.notime <- clx(syrfe.cumvio.notime, 1, alliance$state[!is.na(alliance$lag_num)])
# store model coeffients and clustered se's
coef.matrix[2:6,1] <- clustse.syrfe.cumvio.notime[2:6,1]
se.matrix[2:6,1] <- clustse.syrfe.cumvio.notime[2:6,2]
#variance covariance matrix
vcov.clustse.syrfe.cumvio.notime <- clx2(syrfe.cumvio.notime, 1, alliance$state[!is.na(alliance$lag_num)])
#Model 2
syrfe.rep.notime <- standardize(glm(alliance_gain ~  rep_dum_mem_indef + majpow + democ_dum + lag_num + I(lag_num^2) + as.factor(region1), family= binomial(link=logit),  data=alliance, x=TRUE))
clustse.syrfe.rep.notime <- clx(syrfe.rep.notime, 1, alliance$state[!is.na(alliance$lag_num)])
# store model coeffients and clustered se's
coef.matrix[c(1,3:6),4] <- clustse.syrfe.rep.notime[2:6,1]
se.matrix[c(1,3:6),4] <- clustse.syrfe.rep.notime[2:6,2]
### Perform Clarke test for
### non-nested models on Model's 1 and 2 as described on
### page 29 of the manuscript
clarke(syrfe.cumvio.notime, syrfe.rep.notime)
# Model 3 and 4 comparing our measures with state clustered SE
# plus time dynamics
#Model 3
syrfe.cumvio.time <- standardize(glm(alliance_gain ~  cum_vio_dum+ majpow + democ_dum + lag_num  + I(lag_num^2)  + t + I(t^2) + I(t^3)  + as.factor(region1), family= binomial(link=logit),  data=alliance))
clustse.syrfe.cumvio.time <- clx(syrfe.cumvio.time, 1, alliance$state[!is.na(alliance$lag_num)])
# store model coeffients and clustered se's
coef.matrix[2:6,2] <- clustse.syrfe.cumvio.time[2:6,1]
se.matrix[2:6,2] <- clustse.syrfe.cumvio.time[2:6,2]
#Model 4
syrfe.rep.time <- standardize(glm(alliance_gain ~  rep_dum_mem_indef + majpow + democ_dum + lag_num  + I(lag_num^2)  + t + I(t^2) + I(t^3)  + as.factor(region1), family= binomial(link=logit),  data=alliance))
clustse.syrfe.rep.time <- clx(syrfe.rep.time, 1, alliance$state[!is.na(alliance$lag_num)])
# store model coeffients and clustered se's
coef.matrix[c(1,3:6),5] <- clustse.syrfe.rep.time[2:6,1]
se.matrix[c(1,3:6),5] <- clustse.syrfe.rep.time[2:6,2]
#Model 5 and 6  add a decay term
# Model 5
syrfe.cumvio.decay <- standardize(glm(alliance_gain ~  cum_vio_dum + majpow + democ_dum + lag_num  + I(lag_num^2)  + t + I(t^2) + I(t^3)  + decay_gen + I(decay_gen^2) + I(decay_gen^3)+ as.factor(region1), family= binomial(link=logit),  data=alliance))
clustse.syrfe.cumvio.decay <- clx(syrfe.cumvio.decay, 1, alliance$state[!is.na(alliance$lag_num)])
# store model coeffients and clustered se's
coef.matrix[2:6,3] <- clustse.syrfe.cumvio.decay[2:6,1]
se.matrix[2:6,3] <- clustse.syrfe.cumvio.decay[2:6,2]
# Model 6
# unnormalized model for later
syrfe.rep.decay <- glm(alliance_gain ~  rep_dum_mem_indef + majpow + democ_dum + lag_num  + I(lag_num^2)  + t + I(t^2) + I(t^3)  + decay_gen + I(decay_gen^2) + I(decay_gen^3)+ as.factor(region1), family= binomial(link=logit),  data=alliance)
clustse.syrfe.rep.decay <- clx(syrfe.rep.decay, 1, alliance$state[!is.na(alliance$lag_num)])
# normalized for Figure2
syrfe.rep.decay.norm <- standardize(glm(alliance_gain ~  rep_dum_mem_indef + majpow + democ_dum + lag_num  + I(lag_num^2)  + t + I(t^2) + I(t^3)  + decay_gen + I(decay_gen^2) + I(decay_gen^3)+ as.factor(region1), family= binomial(link=logit),  data=alliance))
clustse.syrfe.rep.decay.norm <- clx(syrfe.rep.decay.norm, 1, alliance$state[!is.na(alliance$lag_num)])
# store model coeffients and clustered se's
coef.matrix[c(1,3:6),6] <- clustse.syrfe.rep.decay.norm[2:6,1]
se.matrix[c(1,3:6),6] <- clustse.syrfe.rep.decay.norm[2:6,2]
# Run additional code to make  Figure 2
source("table_graph_code.R")
######################################
# Create Figure 3 from the main text #
######################################
# This graphs the estimated effect of reputation
attach(alliance)
sims <- 5000
Bhat <- as.vector(clustse.syrfe.rep.decay[,1])
par.draws <- mvrnorm(sims, mu=Bhat, Sigma=clx2(syrfe.rep.decay, 1, alliance$state[!is.na(alliance$lag_num)]))
# set all variables to their median
X1 <- cbind(1, 0, median(majpow, na.rm=T), median(democ_dum, na.rm=T), median(lag_num, na.rm=T), median(lag_num, na.rm=T)^2, median(t, na.rm=T), median(t, na.rm=T)^2, median(t, na.rm=T)^3, median(decay_gen, na.rm=T), median(decay_gen, na.rm=T)^2, median(decay_gen, na.rm=T)^3, 0,0,0,0)
est.fx <- c()
for(i in 6:81){
X1[2] <- i
mu <- 1/(1+exp(-(par.draws%*%t(X1))))
mean.fx <- mean(mu)
sd.fx <- sd(mu)
ci.fx <- quantile(mu, c(.025, 0.975))
est.fx <- rbind(est.fx, c(mean.fx,sd.fx,ci.fx))
}
est.fx <- as.data.frame(est.fx)
names(est.fx) <- c("mean","sd","ci.lo", "ci.hi")
fx.plot <- ggplot(data=est.fx, aes(x=seq(6,81,1), y=mean))+
geom_ribbon(aes(x=seq(6,81,1), ymin=ci.lo, ymax = ci.hi, alpha=.3))+
geom_line(size = 2)+
xlab("Unreliability Score")+ylab("Pr(Alliance Gain | 1 Violation)")
fx.plot <- fx.plot + ggtitle("The Probability of Gaining a new Alliance\nas a Function of Unreliability") +theme_bw() +  theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank(), legend.position="none")
ggsave(fx.plot, file="Figure_3.pdf", width= 6, height =6)
####################################
# Create figure A1 in the appendix #
####################################
timeseq <- 0:max(alliance$t, na.rm=T)
# Set all variables to their median
newX2 <- cbind(1, median(rep_dum_mem_indef, na.rm=T), median(majpow, na.rm=T), median(democ_dum, na.rm=T), median(lag_num, na.rm=T), median(lag_num, na.rm=T)^2, timeseq, timeseq^2,timeseq^3,0,0,0, 0,0,0,0)
Bhat2 <- as.vector(clustse.syrfe.rep.decay[,1])
XBhat2 <- as.matrix(newX2)%*%Bhat2
prob2 <- 1/(1+exp(-XBhat2))
time.plot <- qplot(x=timeseq, y=prob2, geom="line", xlab="t", ylab="Pr(Alliance Gain | t, median(X)")
time.plot <- time.plot + ggtitle("Probability of Gaining a New Alliance\nAfter an Alliance Gain") +theme_bw() +theme( panel.grid.minor = element_blank(), panel.grid.major = element_blank())
ggsave(time.plot, file="Figure_A1.pdf", height=6, width=6)
####################################
# Create Figure A2 in the appendix #
####################################
#graph decay
decayseq <- 0:max(alliance$decay_gen, na.rm=T)
#Set all variables to their median
newX <- cbind(1, median(rep_dum_mem_indef, na.rm=T), median(majpow, na.rm=T), median(democ_dum, na.rm=T), median(lag_num, na.rm=T), median(lag_num, na.rm=T)^2,median(t, na.rm=T), median(t, na.rm=T)^2, median(t, na.rm=T)^3, decayseq, decayseq^2,decayseq^3, 0,0,0,0)
Bhat <- as.vector(clustse.syrfe.rep.decay[,1])
XBhat <- as.matrix(newX)%*%Bhat
prob <- 1/(1+exp(-XBhat))
decay.plot <- qplot(x=decayseq, y=prob, geom="line", xlab="t", ylab="Pr(Alliance Gain | t, median(X)")
decay.plot <- decay.plot + ggtitle("Probability of Gaining a New Alliance\nAfter an Alliance Violation")+theme_bw()+theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank())
ggsave(decay.plot, file="Figure_A2.pdf", width= 6, height=6)
